
C     $Id: bkmp2.F 6 2015-03-03 18:40:49Z wangsl2001@gmail.com $

c
c  _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
c  _/          surface950621 H3 surface evaluation routines          _/
c  _/              (also known as the bkmp2 surface)                 _/
c  _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
c
      subroutine surfid(h3id)
c-----------------------------------------------------------------
c  h3 surface identification number (based on date of fitting)
c-----------------------------------------------------------------
      implicit double precision (a-h,o-z)
c-sngl;	h3id = 950621.0
	h3id = 950621.0d0
      return
      end
c
      subroutine bkmp2(r,Vtot,dVtot,id)
c-----------------------------------------------------------------
c Calculate total H3 potential from all of its parts
c if (id.gt.0) also calculate the dV/dr derivatives
c all distances are in bohrs and all energies are in hartrees
c
c For a discussion of this surface, see:
c    A.I.Boothroyd, W.J.Keogh, P.G.Martin, M.R.Peterson 
c    Journal of Chemical Physics 95 pp4343-4359 (Sept15/91)
c    and a paper in press in JCP (accepted January 1996)
c
c Note: this file contains parameters for a surface refitted
c       on june21/95 to a set of several thousand ab initio
c       points.  the '706' surface parameters have been
c       commented out.  The routine names have been modified
c       slightly (usually a '95' appended) so that a programme
c       could easily call and compare both of our surfaces.
c
c Note: the surface parameter values as published lead to an
c       anomolously deep van der Waals well for a very compact
c       H2 molecule (say r=0.8).  After that paper was submitted,
c       this problem was fixed and the corrected Cbend coefficients
c       are used in this version of the surface (version 706).
c       The old coefficients are still in subr.vbcb but have been
c       commented out.
c
c any QUESTIONS/PROBLEMS/COMMENTS concerning this programme can be
c addressed to :    wkeogh@alchemy.chem.utoronto.ca
c (if necessary, via   boothroy@cita.utoronto.ca    boothroy@utordop.bitnet
c                 or   pgmartin@cita.utoronto.ca    pgmartin@utordop.bitnet )
c
c version:
c apr12/95 ... parameters for surface850308 added
c jul27/91 ... surf706d.out Cbend values put in
c----------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension r(3),dVtot(3)
      dimension dVlon(3),dVas(3),dVbnda(3),dVbndb(3),
     .          dCal(3), dCas(3),dCbnda(3),dCbndb(3)
      dimension vb(2,25),cb(2,25)
c     dimension vbp(2,25),cbp(2,25)    
c     dimension dB1a(3),dB1b(3)
      dimension dT(3),T(3)
c     dimension dB2(3),dB3(6),expass(4)
      ipr = 0

      call chgeom(r,ivalid)
      if(ivalid.lt.1)then
c-sngl;         vtot     = 99.0
         vtot     = 99.0d0
         dvtot(1) = 0.d0
         dvtot(1) = 0.d0
         dvtot(1) = 0.d0
         write(6,6999) vtot,dvtot
         return
      end if
 6999 format('invalid geometry, v,dv=',4f8.2)
c
c
c  zero everything to avoid any 'funny' values:
      Vlon  = 0.d0
      Vas   = 0.d0
      Vbnda = 0.d0
      Vbndb = 0.d0
      Cal   = 0.d0
      Cas   = 0.d0
      cbnda = 0.d0
      cbndb = 0.d0
      do i=1,3
	 dVlon(i) = 0.d0
	 dVas(i)  = 0.d0
	 dVbnda(i)= 0.d0
	 dVbndb(i)= 0.d0
	 dCal(i)  = 0.d0
	 dCas(i)  = 0.d0
	 dCbnda(i)= 0.d0
	 dCbndb(i)= 0.d0
      end do  
c  zero the vb and cb arrays (used only for numerical derivatives)
      do i=1,25
	 vb(1,i) = 0.d0
	 vb(2,i) = 0.d0
	 cb(1,i) = 0.d0
	 cb(2,i) = 0.d0
      end do  
      call h3lond95( r, Vlon, dVlon )
      call vascal95( r, Vas, dVas )
c    now do any corrections required for compact geometries:
      call compac95(r,icompc,T,dT)
c    oct.3/90 compact routines only called for compact geometries:
      if(icompc.ge.1)then
         call csym95  ( r, cal, dCal )
         call casym95 ( r, Cas, dCas,  1 ,T,dT)
      end if
      call vbcb95(r,icompc,T,dT,ipr,
     .          Vbnda,Vbndb,dVbnda,dVbndb,
     .          Cbnda,Cbndb,dCbnda,dCbndb)
c    add up the various parts of the potential:
      Vtot = Vlon + Vas + Vbnda + Vbndb + Cal + Cas + Cbnda + Cbndb
c    add up the various parts of the derivative:
      dVtot(1) =   dVlon(1) + dVas(1) + dVbnda(1) + dVbndb(1) 
     .           +  dCal(1) + dCas(1) + dCbnda(1) + dCbndb(1)
      dVtot(2) =   dVlon(2) + dVas(2) + dVbnda(2) + dVbndb(2) 
     .           +  dCal(2) + dCas(2) + dCbnda(2) + dCbndb(2)
      dVtot(3) =   dVlon(3) + dVas(3) + dVbnda(3) + dVbndb(3) 
     .           +  dCal(3) + dCas(3) + dCbnda(3) + dCbndb(3)

      if(ipr.gt.0) then
	  write(7,*)
	  write(7,*) 'h3tot enter --------------------------------'
          write(7,7000) r,icompc,vtot,
     .             vlon,vas,cal,cas,vbnda,vbndb,cbnda,cbndb
	  if(ipr.gt.1)then
		   write(7,7100) dVtot,dVlon,dVas,dVbnda,dVbndb,
     .                           dCal,dCas,dCbnda,dCbndb
	  end if
	  write(7,7400) (vb(1,i),i=1,5),(vb(2,i),i=1,5)
	  write(7,7500) (cb(1,i),i=1,9),(cb(2,i),i=1,9)
          write(7,7999) Vbnda,Vbndb,Cbnda,Cbndb,  
     .                  dVbnda,dvbndb,dCbnda,dCbndb
	  write(7,*) 'exiting subr.h3tot'
	  write(7,*) 'h3tot exit ---------------------------------'
      end if
      return
 7400 format('Vba values: ',5(1x,g12.6),/,'Vbb values: ',5(1x,g12.6))
 7500 format('Cba values: ',3(1x,f16.8),/,
     .       '            ',3(1x,f16.8),/,
     .       '            ',3(1x,f16.8),/,
     .       'Cbb values: ',3(1x,f16.8),/,
     .       '            ',3(1x,f16.8),/,
     .       '            ',3(1x,f16.8))  
 7000 format(5x,'    r =',3(1x,f16.10),' icompac = ',i1,/,
     .       5x,'Vtot  = ',f18.12,/,
     .       5x,'Vlon  = ',f18.12,'       Vas   = ',g18.12,/,
     .       5x,'Cal   = ',g18.12,'       Cas   = ',g18.12,/,
     .       5x,'Vbnda = ',g18.12,'       Vbndb = ',g18.12,/,
     .       5x,'Cbnda = ',g18.12,'       Cbndb = ',g18.12)
 7100 format(' dVtot   = ',3(1x,g18.12),/,
     .       ' dVlon   = ',3(1x,g18.12),/,
     .       ' dVas    = ',3(1x,g18.12),/,
     .       ' dVbnda  = ',3(1x,g18.12),/,
     .       ' dVbndb  = ',3(1x,g18.12),/,
     .       ' dCal    = ',3(1x,g18.12),/,
     .       ' dCas    = ',3(1x,g18.12),/,
     .       ' dCbnda  = ',3(1x,g18.12),/,
     .       ' dCbndb  = ',3(1x,g18.12))
 7750 format(3(1x,g16.10))
 7751 format(/,3(1x,g16.10))
 7999 format('   Vbnda        Vbndb        Cbnda        Cbndb ',/,
     .        4(e12.6,2x),/,
     .       'dVbnda: ',3(e16.10,1x),/,
     .       'dVbndb: ',3(e16.10,1x),/,
     .       'dCbnda: ',3(e16.10,1x),/,
     .       'dCbndb: ',3(e16.10,1x))
      end
c
      subroutine triplet95(r,E3,id)
C--------------------------------------------------------------------
c apr12/95 surface950308 values added
c oct04/90 surface626 values added
C  H2 triplet curve and derivatives:
c  calculates triplet potential and first derivative
c  uses truhlar horowitz equation with our extension
c  uses the Johnson correction at short distances (r < rr)
c     if r .ge. rr         use modified t/h triplet equation
c     if r .le. rl         use the Johnson correction
c     in between           use the transition equation
C--------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension E3(3),E1(3)
      parameter( rl=0.95d0, rr=1.15d0 )
      parameter( z1=1.d0, z2=2.d0, z4=4.d0)
c triplet and johnson values mar08/95                                                       
c     parameter( 
c    .  a1= -0.0460730469,a2=-16.3948753611,a3=-27.6026451962,
c    .  a4=  2.0383359014,a5= -6.8887433326,a6=  1.6216005938,
c    .  C1= -0.4110341110310669,C2= -0.0767398928636804,
c    .  C3=  0.4302169435606523)
c triplet values jun21                                                                
c-sngl;      parameter( 
c-sngl;     .  a1= -0.0298546962,a2=-23.9604445036,a3=-42.5185569474,
c-sngl;     .  a4=  2.0382390988,a5=-11.5214861455,a6=  1.5309487826,
c-sngl;     .  C1= -0.4106358351531854,C2= -0.0770355790707090,
c-sngl;     .  C3=  0.4303193846943223) !jun21 fit
      parameter( 
     .  a1= -0.0298546962d0,a2=-23.9604445036d0,a3=-42.5185569474d0,
     .  a4=  2.0382390988d0,a5=-11.5214861455d0,a6=  1.5309487826d0,
     .  C1= -0.4106358351531854d0,C2= -0.0770355790707090d0,
     .  C3=  0.4303193846943223d0) !jun21 fit
c  triplet and johnson values from fit621:
c     parameter(a1=-0.0253496194,a2=-29.2302126444,
c    .          a3=-50.7225015503,a4= 2.0452676876 ,
c    .          a5=-12.2408908509,a6= 1.6733157383 )
c     parameter( C1=-0.4170298146519658,
c    .           C2=-0.0746027774843370,C3= 0.4297899952237434 )
c surface parameters from fit601.out
c     parameter(a1=-0.66129429,a2=-1.99434198,a3=-2.37604328,
c    .          a4= 2.08107802,a5=-0.0313032510,a6=3.76546699,
c    .          C1=-0.4222590135447196,C2=-0.0731117796738824,
c    .          C3= 0.4295918082189010 )
      ipr = 0
      E3(3) = 0.d0
      if(r.ge.rr )then
c       modified truhlar/horowitz triplet equation:
         exdr = dexp( -a4*r )
         rsq  = r*r
         ra6  = r**(-a6)
         E3(1) = a1* ( a2 + r + a3*rsq + a5*ra6 )*exdr 
c       first derivative of triplet curve:
         ra61 = r**(-a6-z1)
         E3(2) = a1*exdr*
     .   ( z1 -a2*a4 +(z2*a3-a4)*r -a3*a4*rsq -a5*a6*ra61 -a4*a5*ra6 )
      end if
      if(r.lt.rr )then
          dr = r- rl  
          call vh2opt95(r,e1,2)
	  if( r.le.rl ) then
c          Johnson triplet equation:
            E3(1) = E1(1) + c2*dr + c3
	    E3(2) = E1(2) + c2
          else                 
c          Transition equation:    
            E3(1) = E1(1) + c1*dr*dr*dr + c2*dr + c3
	    E3(2) = E1(2) + 3.d0*c1*dr*dr + c2
	  end if
      end if
      if(ipr.gt.0) write(7,7100) e3
 7100 format('  triplet:  E3 = ',3(1x,f12.8))
      return
      end
C
c
      subroutine vH2opt95(r,E,ideriv)  
C-----------------------------------------------------------------
c jul09/90 ... super duper speedy version
C self-contained version of schwenke's H2 potential
C all distances in bohrs and all energies in Hartrees
C (1st deriv added on May 2 1989; 2nd deriv on May 28 1989)
C-----------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension E(3)
c-sngl;      parameter(a0=0.03537359271649620,    a1= 2.013977588700072    ,
c-sngl;     .  a2= -2.827452449964767        ,    a3= 2.713257715593500    ,
c-sngl;     .  a4= -2.792039234205731        ,    a5= 2.166542078766724    ,
c-sngl;     .  a6= -1.272679684173909        ,    a7= 0.5630423099212294   ,
c-sngl;     .  a8= -0.1879397372273814       ,    a9= 0.04719891893374140  ,
c-sngl;     . a10= -0.008851622656489644     ,   a11= 0.001224998776243630 ,
      parameter(
     .  a0=0.03537359271649620d0    ,    a1= 2.013977588700072d0    ,
     .  a2= -2.827452449964767d0    ,    a3= 2.713257715593500d0    ,
     .  a4= -2.792039234205731d0    ,    a5= 2.166542078766724d0    ,
     .  a6= -1.272679684173909d0    ,    a7= 0.5630423099212294d0   ,
     .  a8= -0.1879397372273814d0   ,    a9= 0.04719891893374140d0  ,
     . a10= -0.008851622656489644d0 ,   a11= 0.001224998776243630d0 ,
     . a12= -1.227820520228028d-04  ,   a13= 8.638783190083473d-06  ,
     . a14= -4.036967926499151d-07  ,   a15= 1.123286608335365d-08  ,
     . a16= -1.406619156782167d-10 )
      parameter( r0=3.5284882d0,  DD=0.160979391d0,
     .           c6=6.499027d0,   c8=124.3991d0,    c10=3285.828d0)
c     Ediss = 0.174445d0
      E(1)  = -999.d0
      E(2)  = -999.d0
      E(3)  = -999.d0
      r2  = r  *r
      r3  = r2 *r
      r4  = r3 *r
      r5  = r4 *r
      r6  = r5 *r
      r7  = r6 *r
      r8  = r7 *r
      r9  = r8 *r
      r10 = r9 *r
      r11 = r10*r
      r12 = r11*r
      r13 = r12*r
      r14 = r13*r
      r15 = r14*r
      r02 = r0*r0
      r04 = r02*r02
      r06 = r04*r02
      rr2  = r2 + r02
      rr4  = r4 + r04
      rr6  = r6 + r06
      rr25 = rr2*rr2*rr2*rr2*rr2
c     general term:  a(i)*r(i-1),  i=0,16
      alphaR =  a0/r + a1
     .            + a2 *r   + a3 *r2  + a4 *r3  + a5 *r4  + a6 *r5
     .            + a7 *r6  + a8 *r7  + a9 *r8  + a10*r9  + a11*r10
     .            + a12*r11 + a13*r12 + a14*r13 + a15*r14 + a16*r15
      exalph = dexp(alphaR)
      vsr = DD*(exalph-1.d0)*(exalph-1.d0) - DD
      vlr =  -C6/rr6   -C8/(rr4*rr4)   -C10/rr25
      E(1)=  vsr +  vlr 
c
C  calculate first derivative if required:
      if(ideriv.ge.1)then        
         r3 = r2*r
         r5 = r4*r
         rr26 = rr25*rr2
         rr43 = rr4*rr4*rr4
c        general term:  (i-1)*a(i)*r**(i-2)  ,  i=0,16
	 dalphR = -a0/r2 + a2 + 2.d0*a3*r
     .          +  3.d0 *a4 *r2  +  4.d0 *a5 *r3  +  5.d0 *a6 *r4 
     .          +  6.d0 *a7 *r5  +  7.d0 *a8 *r6  +  8.d0 *a9 *r7 
     .          +  9.d0 *a10*r8  + 10.d0 *a11*r9  + 11.d0 *a12*r10
     .          + 12.d0 *a13*r11 + 13.d0 *a14*r12 + 14.d0 *a15*r13
     .          + 15.d0 *a16*r14
	 dvsr = 2.d0 *DD *(exalph-1.d0) *exalph *dalphR
	 dvlr =    6.d0*C6*r5 / (rr6*rr6)
     .          +  8.d0*C8*r3 / rr43  + 10.d0*C10*r / rr26
	 E(2) = dvsr + dvlr
      end if
c
C  calculate second derivative if required:
      if(ideriv.ge.2)then
         r10  = r6*r4
         rr27 = rr26*rr2
         rr44 = rr43*rr4
         rr62 = rr6*rr6
         rr63 = rr62*rr6
c        general term: (i-1)*(i-2)*a(i)*r**(i-3),    i=0,16
	 ddalph = 2.d0*a0/r3 +  2.d0*a3      +  6.d0*a4 *r
     .     + 12.d0*a5 *r2    + 20.d0*a6 *r3  + 30.d0*a7 *r4
     .     + 42.d0*a8 *r5    + 56.d0*a9 *r6  + 72.d0*a10*r7
     .     + 90.d0*a11*r8    +110.d0*a12*r9  +132.d0*a13*r10
     .     +156.d0*a14*r11   +182.d0*a15*r12 +210.d0*a16*r13
	 ddvsr = 2.d0 *DD *exalph 
     .   *( (2.d0*exalph-1.d0)*dalphR*dalphR + (exalph-1.d0)*ddalph )
	 ddvlr =- 72.d0* C6 *r10 / rr63 -96.d0 *C8 *r6 / rr44  
     .		-120.d0*C10 *r2  / rr27 +30.d0 *C6 *r4 / rr62
     .		+ 24.d0* C8 *r2  / rr43 +10.d0 *C10    / rr26
	 E(3) = ddvsr + ddvlr
      end if
      return
      end
c
      subroutine h3lond95(r,Vlon,dVlon)
c----------------------------------------------------------------------
c version of may 12/90  ... derivatives corrected (0.5 changed to 0.25)
c calculates the h3 london terms and derivatives 
c modified oct 7/89 to include eps**2 term which rounds off the
c cusp in the h3 potential which occurs at equilateral triangle 
c configurations
c----------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      real*8 Q(3),J(3),Jt
      dimension r(3),e1(3),e3(3),esing(3),etrip(3),dVlon(3)
      dimension de1(3),de3(3)
c     dimension dJ1(3),dJ2(3),dJ3(3)
      parameter( half=0.5d0, two=2.d0 , eps2=1.d-12 )
      ipr = 0
      do i=1,3
         call vh2opt95(r(i),esing,2)
	  e1(i) = esing(1)
	 de1(i) = esing(2)
	 call triplet95(r(i),etrip,2)
	 if(ipr.gt.0) write(7,7400) i,esing,etrip
	  e3(i) = etrip(1)
	 de3(i) = etrip(2)
         q(i)  = half*(E1(i) + E3(i))  
         j(i)  = half*(E1(i) - E3(i)) 
      end do   
      sumQ  =   q(1) + q(2) + q(3)
      sumJ  =   dabs( j(2)-j(1) )**2 
     .        + dabs( j(3)-j(2) )**2 
     .        + dabs( j(3)-j(1) )**2 
      Jt     = half*sumJ + eps2
      rootJt = dsqrt(Jt)
      Vlon   = sumQ - rootJt   
      if(ipr.gt.0) then
	 write(7,7410) sumq,sumj
         write(7,7420) vlon,rootJt
      end if
c  calculate the derivatives with respect to r(i):
      dVlon(1) = half*(dE1(1)+de3(1))
     .         - 0.25d0*(two*j(1)-j(2)-j(3))*(dE1(1)-de3(1))/rootJt
      dVlon(2) = half*(dE1(2)+de3(2))
     .         - 0.25d0*(two*j(2)-j(3)-j(1))*(dE1(2)-de3(2))/rootJt
      dVlon(3) = half*(dE1(3)+de3(3))
     .         - 0.25d0*(two*j(3)-j(1)-j(2))*(dE1(3)-de3(3))/rootJt
      if(ipr.gt.0) then
         write(7,7000) r,e1,e3,vlon
	 write(7,7100) q,j
	 write(7,7200) dVlon
      end if
 7000 format('             r = ',3(1x,f12.6),/,
     .       '            E1 = ',3(1x,f12.8),/,
     .       '            E3 = ',3(1x,f12.8),/,
     .       '          vlon = ',1x,f12.8)
 7100 format(13x,'q = ',3(1x,e12.6),/,13x,'j = ',3(1x,e12.6))
 7200 format('         dVlon = ',3(1x,g12.6))
 7400 format('from subr.london: ',/,
     .       '  using r',i1,':   Esinglet=',3(1x,f12.8),/,
     .       '         ',1x,'    Etriplet=',3(1x,f12.8))
 7410 format('         sumQ = ',g12.6,'        sumJ = ',g12.6)
 7420 format('         Vlon = ',f12.8,'      rootJt = ',g12.6)
      return
      end
c
      subroutine Vascal95(rpass,Vas,dVas)
c------------------------------------------------------------------
c version of apr12/95  950308 values      
c version of oct11/90  fit632c.out values
c version of oct5/90   surf626 values
c  calculate the asymmetric correction term and its derivatives
c  see equations [14] to [16] of truhlar/horowitz 1978 paper
c------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension dVas(3),dA(3),dS(3),rpass(3)
c vasym values jun21                                                            
c-sngl;      parameter( 
c-sngl;     . aa1=0.3788951192E-02, aa2=0.1478100901E-02,
c-sngl;     . aa3=-.1848513849E-03, aa4=0.9230803609E-05,
c-sngl;     . aa5=-.1293180255E-06, aa6=0.5237179303E+00,
c-sngl;     . aa7=-.1112326215E-02) !jun21 fit
      parameter( 
     . aa1=0.3788951192d-02, aa2=0.1478100901d-02,
     . aa3=-.1848513849d-03, aa4=0.9230803609d-05,
     . aa5=-.1293180255d-06, aa6=0.5237179303d+00,
     . aa7=-.1112326215d-02) !jun21 fit
c--- vasym values mar08.95                                                         
c     parameter( 
c    . aa1=0.3759731624E-02, aa2=0.1476254095E-02,
c    . aa3=-.1866759453E-03, aa4=0.9218646237E-05,
c    . aa5=-.1287906069E-06, aa6=0.5201790843E+00,
c    . aa7=-.1062909514E-02)
c--- Vasym values from fit632c ----------------------
c     parameter(
c    . aa1=0.3438222224E-02,aa2=0.1398145763E-02,
c    . aa3=-.1923999449E-03,aa4=0.9712737075E-05,
c    . aa5=-.1263794562E-06,aa6=0.5181432712E+00,
c    . aa7=-.9487002995E-03 )
      ipr = 0
      r1 = rpass(1)
      r2 = rpass(2)
      r3 = rpass(3)
      R  = r1 + r2 + r3
      Rsq = R*R
      Rcu = Rsq*R
C   calculate the Vas term first (eq.14 of truhlar/horowitz)
      call acalc95(r1,r2,r3,A,dA)
        A2 = A *A
        A3 = A2*A
        A4 = A3*A
        A5 = A4*A
        exp1 = dexp(-aa1*Rcu)
	exp6 = dexp(-aa6*R)
        S = aa2*A2 + aa3*A3 + aa4*A4 + aa5*A5
      Vas = S*exp1 +  aa7*A2 *exp6 / R
         dS(1)  = ( 2.d0*aa2*A  +3.d0*aa3*A2 
     .             +4.d0*aa4*A3 +5.d0*aa5*A4) * dA(1)
         dVas(1) =  -3.d0*aa1*Rsq*S*exp1 + dS(1)*exp1
     .             -aa7*A2*exp6/Rsq  + 2.d0*aa7*A*dA(1)*exp6/R
     .             -aa6*aa7*A2*exp6/R
         dS(2)  = ( 2.d0*aa2*A  +3.d0*aa3*A2 
     .             +4.d0*aa4*A3 +5.d0*aa5*A4) * dA(2)
         dVas(2) =  -3.d0*aa1*Rsq*S*exp1 + dS(2)*exp1
     .             -aa7*A2*exp6/Rsq  + 2.d0*aa7*A*dA(2)*exp6/R
     .             -aa6*aa7*A2*exp6/R
         dS(3)  = ( 2.d0*aa2*A  +3.d0*aa3*A2 
     .             +4.d0*aa4*A3 +5.d0*aa5*A4) * dA(3)
         dVas(3) =  -3.d0*aa1*Rsq*S*exp1 + dS(3)*exp1
     .             -aa7*A2*exp6/Rsq  + 2.d0*aa7*A*dA(3)*exp6/R
     .             -aa6*aa7*A2*exp6/R
      if(ipr.gt.0) write(7,7100) Vas,dS,dVas
 7100 format('    from subr.vascal:   Vas  = ',  1x,g12.6, /,
     .       '                        dS   = ',3(1x,g12.6),/,
     .       '                        dVas = ',3(1x,g12.6))
      return
      end
c
      subroutine acalc95(r1,r2,r3,A,dA)
c---------------------------------------------------------------------
c  version of may 1, 1990
c  based on equations from notes date april 4, 1990
c  calculate the A value and its derivatives
c           A  = dabs[ (r1-r2)*(r2-r3)*(r3-r1) ]
c---------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension dA(3)
      ipr = 0
      A = (r1-r2)*(r2-r3)*(r3-r1)
      dA(1) = ( -2.d0*r1 + r2 + r3 )*(r2-r3)
      dA(2) = ( -2.d0*r2 + r3 + r1 )*(r3-r1)
      dA(3) = ( -2.d0*r3 + r1 + r2 )*(r1-r2)
      if(A.lt.0.d0)then
         A = -A
         dA(1) = -dA(1)
         dA(2) = -dA(2)
         dA(3) = -dA(3)
      end if
      if(ipr.gt.0) write(7,7100) A,dA
 7100 format('Acalc enter ---------------------------------------',/,
     .       '     A = ',f15.10,'   dA = ',3(1x,g12.6),/,
     .       'Acalc exit ----------------------------------------')
      return
      end 
c
      subroutine compac95(r,icompc,T,dT)
C---------------------------------------------------------------
c  version of may 14/90 ... calculates T and dT values also
c  decide whether or not this particular geometry is compact,
c  that is, are any of the three distances smaller than the
c  rr value from the johnson correction.
c-----------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension r(3),T(3),dT(3)
      parameter( rr = 1.15d0, rp = 1.25d0 )
c  first see if this is a compact geometry:
      do i=1,3
	 T(i) = 0.d0
	 dT(i)= 0.d0
      end do   
      icompc = 0
      if(r(1).lt.rr) icompc = icompc + 1
      if(r(2).lt.rr) icompc = icompc + 1
      if(r(3).lt.rr) icompc = icompc + 1
      if(icompc.eq.0) return
c calculate the T(i) values:
      do i=1,3
      if(r(i).lt.rr) then
	 top = rr-r(i)
	 bot = rp-r(i)
	 top2 = top  * top
	 top3 = top2 * top
	 bot2 = bot  * bot
         T(i)  = top3/bot
         dT(i) = -3.d0*top2/bot + top3/bot2
      end if
      end do   
      return
      end
C
      subroutine casym95( r,Cas,dCas,id,T,dT )
c---------------------------------------------------------------
c apr12/95 ... surface950308 parameters added
c version of sept14/90 
c  the compact asymmetric correction term and derivatives
c---------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension r(3),T(3)
      dimension dPr(3),dT(3),dsumT(3),dterm1(3),determ(3),dseries(3),
     .          dCas(3),dA(3)
c Casym values jun21                                                            
c-sngl;      parameter( 
c-sngl;     . u1=0.2210243144E+00, u2=0.4367417579E+00, u3=0.6994985432E-02,
c-sngl;     . u4=0.1491096501E+01, u5=0.1602896673E+01, u6=-.2821747323E+01,
c-sngl;     . u7=0.4948310833E+00, u8=-.3540394679E-01, u9=-.3305809954E+01,
c-sngl;     .u10=0.3644382172E+01,u11=-.9997570970E+00,u12=0.7989919534E-01,
c-sngl;     .u13=-.1075807322E-02) !jun21 fit
      parameter( 
     . u1=0.2210243144d+00, u2=0.4367417579d+00, u3=0.6994985432d-02,
     . u4=0.1491096501d+01, u5=0.1602896673d+01, u6=-.2821747323d+01,
     . u7=0.4948310833d+00, u8=-.3540394679d-01, u9=-.3305809954d+01,
     .u10=0.3644382172d+01,u11=-.9997570970d+00,u12=0.7989919534d-01,
     .u13=-.1075807322d-02) !jun21 fit
c---- Casym values mar08.95                                                         
c     parameter( 
c    . u1=0.5120831287E+00, u2=0.1002277242E+01, u3=0.6850007419E-02,
c    . u4=-.2038751706E+01, u5=0.7027811909E+01, u6=-.4881767278E+01,
c    . u7=0.8801769106E+00, u8=-.6296419648E-01, u9=-.8125516783E+01,
c    .u10=0.6073964424E+01,u11=-.1451402523E+01,u12=0.1165084183E+00,
c    .u13=-.1176579871E-02)
c---- Casym values from fit633a   
c     parameter( 
c    . u1=0.1537481166E+00, u2=0.2745950036E+00, u3=0.7501206780E-02,
c    . u4=0.3119136023E+01, u5=0.9969170798E+00, u6=-.3373682823E+01,
c    . u7=0.6807215913E+00, u8=-.4920325491E-01, u9=-.3919467989E+01,
c    .u10=0.5085532326E+01,u11=-.1415264778E+01,u12=0.1138681785E+00,
c    .u13=-.1525367566E-02)
      ipr = 0
      cas = 0.d0
      dCas(1) = 0.d0
      dCas(2) = 0.d0
      dCas(3) = 0.d0
      call acalc95(r(1),r(2),r(3),A,dA)
      A2 = A*A
c     if(A.eq.0.d0) return
      sumT = T(1) + T(2) + T(3)
      Sr   = r(1) + r(2) + r(3)
      Pr   = r(1) * r(2) * r(3)
      Sr2  = Sr*Sr
      Sr3  = Sr2*Sr
      Pr2  = Pr*Pr
      Pr3  = Pr2*Pr
c    write out the series explicitly:
      series = 1.d0 + u4/Pr2 +  u5/Pr + u6 +   u7*Pr +  u8*Pr2
     .           + A*(u9/Pr2 + u10/Pr + u11 + u12*Pr + u13*Pr2)
      term1  = u1/Pr**u2
      eterm  = dexp(-u3*Sr3)
      Cas    = sumt * a2 * term1 * series * eterm
      if(ipr.gt.0) then
	    write(7,*)
            write(7,*) 'Casym: ------------------'
	    write(7,*) ' T(i) = ',T
	    write(7,*) 'dT(i) = ',dT
	    write(7,*) ' id = ',id
	    write(7,7400) series,term1,eterm,Cas
	    write(7,*) ' Cas = ',Cas
      end if
      if(id.gt.0)then 
	    dPr(1) = r(2)*r(3)
	    dPr(2) = r(3)*r(1)
	    dPr(3) = r(1)*r(2)
	 do i=1,3
	    dsumt(i) = dt(i)
	    dterm1(i)=  -1.d0*u1*u2*Pr**(-u2-1.d0)*dPr(i)
	    determ(i)= eterm *(-3.d0*u3*Sr2)
	    dseries(i)=
     .        dPr(i)*( -2.d0*u4/Pr3  -u5/Pr2 +u7  +2.d0*u8*Pr    )
     .        +dA(i)*( u9/Pr2 +u10/Pr +u11 +u12*Pr +u13*Pr2    )
     .        +A*dPr(i)*( -2.d0*u9/Pr3 -u10/Pr2 +u12 + 2.d0*u13*Pr )
	    dCas(i) = 
     .        dsumt(i)      * a2    * term1 * series * eterm
     .      + 2.d0*a*dA(i)  * sumt  * term1 * series * eterm
     .      + dterm1(i)     * sumt  * a2    * series * eterm
     .      + dseries(i)    * sumt  * a2    * term1  * eterm
     .      + determ(i)     * sumt  * a2    * term1  * series 
         end do   
	 if(ipr.gt.0) then
	    write(7,7500) dPr,dt,dterm1,determ,dseries,dCas
	    write(7,*) ' dCas = ',dCas
	 end if
      end if
      return
 7400 format('subr.Cas:   series = ',g12.6,'     term1 = ',g12.6,/,
     .       '             eterm = ',g12.6,'      Cas  = ',g12.6)
 7500 format('    derivatives: dPr = ',3(1x,g12.6),/,
     .       '                 dT  = ',3(1x,g12.6),/,
     .       '              dterm1 = ',3(1x,g12.6),/,
     .       '              determ = ',3(1x,g12.6),/,
     .       '             dseries = ',3(1x,g12.6),/,
     .       '              dCas   = ',3(1x,g12.6))  
      end
c
      subroutine csym95(r,Cal,dCal)
c---------------------------------------------------------------
c version of apr12/95  950308 values 
c version of oct12/90  surf636 values
c  calculate the 'compact all' correction term and derivatives
c  a correction term (added sept 11/89), which adds a small
c  correction to all compact geometries
c---------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension r(3),dCal(3),G(3)
c     dimension t(3)
      dimension dG(3),sumv(3)
      parameter( rr=1.15d0, rp=1.25d0 )
c Csym values jun21                                                             
c-sngl;      parameter( 
c-sngl;     . v1=-.2071708868E+00, v2=-.5672350377E+00, v3=0.9058780367E-02) !jun21 fit
      parameter( 
     . v1=-.2071708868d+00, v2=-.5672350377d+00, v3=0.9058780367d-02) !jun21 fit
c---- Csym values mar08/95                                                          
c     parameter( 
c    . v1=-.2210049400E+00, v2=-.7054469608E+00, v3=0.4088898394E-02)
c---- Csym values from fit633a                                                  
c     parameter( 
c    . v1=-.2070776049E+00, v2=-.5350898737E+00, v3=0.1011861942E-01)
      cal   = 0.d0
      ipr   = 0
      Sr    = r(1)+r(2)+r(3)
      Sr2   = Sr*Sr
      Sr3   = Sr*Sr2  
      exp3  = dexp( -v3*Sr3 )
      dexp3 = -3.d0*v3*Sr2*exp3
      do i=1,3
         ri      = r(i)
	 rrri    = rr-ri
	 rrri2   = rrri*rrri
	 rrri3   = rrri*rrri2
	 rpri    = rp-ri
	 rpri2   = rpri*rpri
         G(i)    = 0.d0
	 dG(i)   = 0.d0
	 sumv(i) = v1+v1*v2*ri
         if(ri.lt.rr) then
	    G(i)  = (rrri3/rpri) *sumv(i)
  	    dG(i) = (rrri3/rpri2)*sumv(i) 
     .                -3.d0*(rrri2/rpri)*sumv(i)
     .                + (rrri3/rpri)*v1*v2
	 end if
      end do   
      sumG = G(1) + G(2) + G(3)
      cal  = sumG*exp3
      dCal(1) = dG(1)*exp3 + sumG*dexp3
      dCal(2) = dG(2)*exp3 + sumG*dexp3
      dCal(3) = dG(3)*exp3 + sumG*dexp3
      if(ipr.gt.0)then
	 write(7,*)
         write(7,7100) Cal
	 write(7,7000) G,dG,sumv,dCal
	 write(7,*) 'exiting subr.csym'
      end if
 7100 format('Csym enter ----------------------------------------',/,
     .       '    Cal = ',1x,g20.14)
 7000 format('      G = ',3(1x,g12.6),/,
     .       '     dG = ',3(1x,g12.6),/,
     .       '   sumv = ',3(1x,g12.6),/,
     .       '   dCal = ',3(1x,g16.10),/,
     .       'Csym exit -----------------------------------------')
      return
      end
c
      subroutine vbcb95(rpass,icompc,T,dT,ipr,
     .                  Vbnda,Vbndb,dVbnda,dVbndb,
     .                  Cbnda,Cbndb,dCbnda,dCbndb)
C---------------------------------------------------------------
c apr12/95 surface950308 values added
c jul24/91 fit705 cbend values added
c feb27/91 modified to match equation in h3 paper more closely
c nov 4/90 Vbend coefficients now a,g  Cbend coeff's still c,d
c in this version, the derivatives are always calculated
c subroutines vbend, cbend and B1ab all combined into this one
c module in order to improve efficiency by not having to pass
c around B1a,B1b,B2,B3a,B3b functions and derivatives
c  B1a = 1 - sum of [ P1(cos(theta(i))) ]
c  B1b = 1 - sum of [ P3(cos(theta(i))) ]
c--------------------------------------------------------------
      implicit double precision (a-h,o-z)
c     arrays required for B1,B2,B3 calculations:
      dimension rpass(3),dB1a(3),dB1b(3),dB2(3),dB3(6)
c     dimension th(3)
c     arrays required for Vbend calculations:
      dimension dVbnda(3),dVbndb(3)
      dimension dVb1(3),dvb2(3),dvb3(3),dvb4(3),dvb5(3)
c     dimension dVbnd(3)
      dimension Vbnd(2)
      dimension vb(2,25)   
c     arrays required for Cbend calculations:
      dimension Cbnds(2),dCbnds(2,3)
      dimension dCb1(3),dCb2(3),dCb3(3),dCb4(3),dCb5(3),dCb6(3),
     .          dCb7(3),dCb8(3)
      dimension T(3),dP(3)
      dimension dT(3),dCbnda(3),dCbndb(3)
      dimension cb(2,25)   
      parameter( z58=0.625d0, z38=0.375d0 )  
c     parameters required for Vbend calculations:
      parameter(beta1=0.52d0, beta2=0.052d0, beta3=0.79d0 ) 
c vbenda terms jun21         expterms= 0.52d0, 0.052d0, 0.79d0                        
c-sngl;      parameter( 
c-sngl;     .a11=-.1838073394E+03,a12=0.1334593242E+02,a13=-.2358129537E+00,
c-sngl;     .a21=-.4668193478E+01,a22=0.7197506670E+01,a23=0.2162004275E+02,
c-sngl;     .a24=0.2106294028E+02,a31=0.4242962586E+01,a32=0.4453505045E+01,
c-sngl;     .a41=-.1456918088E+00,a42=-.1692657366E-01,a43=0.1279520698E+01,
c-sngl;     .a44=-.4898940075E+00,a51=0.1742295219E+03,a52=0.3142175348E+02,
c-sngl;     .a53=0.5152903406E+01) !jun21 fit
      parameter( 
     .a11=-.1838073394d+03,a12=0.1334593242d+02,a13=-.2358129537d+00,
     .a21=-.4668193478d+01,a22=0.7197506670d+01,a23=0.2162004275d+02,
     .a24=0.2106294028d+02,a31=0.4242962586d+01,a32=0.4453505045d+01,
     .a41=-.1456918088d+00,a42=-.1692657366d-01,a43=0.1279520698d+01,
     .a44=-.4898940075d+00,a51=0.1742295219d+03,a52=0.3142175348d+02,
     .a53=0.5152903406d+01) !jun21 fit
c vbendb terms jun21                                                                 
c-sngl;      parameter( 
c-sngl;     .g11=-.4765732725E+02,g12=0.3648933563E+01,g13=-.7141145244E-01,
c-sngl;     .g21=0.1002349176E-01,g22=0.9989856329E-02,g23=-.4161953634E-02,
c-sngl;     .g24=0.9075807910E-03,g31=-.2693628729E+00,g32=-.1399065763E-01,
c-sngl;     .g41=-.1417634346E-01,g42=-.4870024792E-03,g43=0.1312231847E+00,
c-sngl;     .g44=-.4409850519E-01,g51=0.5382970863E+02,g52=0.4587102824E+01,
c-sngl;     .g53=0.1768550515E+01) !jun21 fit
      parameter( 
     .g11=-.4765732725d+02,g12=0.3648933563d+01,g13=-.7141145244d-01,
     .g21=0.1002349176d-01,g22=0.9989856329d-02,g23=-.4161953634d-02,
     .g24=0.9075807910d-03,g31=-.2693628729d+00,g32=-.1399065763d-01,
     .g41=-.1417634346d-01,g42=-.4870024792d-03,g43=0.1312231847d+00,
     .g44=-.4409850519d-01,g51=0.5382970863d+02,g52=0.4587102824d+01,
     .g53=0.1768550515d+01) !jun21 fit
c    vbenda mar08/95   expterms= 0.52d0, 0.052d0, 0.79d0
c     parameter( 
c    .a11=-.1779469255E+03,a12=0.1292511538E+02,a13=-.2284450520E+00,
c    .a21=-.4671094689E+01,a22=0.7810423901E+01,a23=0.2359968959E+02,
c    .a24=0.2293876132E+02,a31=0.4191795902E+01,a32=0.4456025797E+01,
c    .a41=-.1419862958E+00,a42=-.1667760578E-01,a43=0.1246997125E+01,
c    .a44=-.4777517290E+00,a51=0.1681125268E+03,a52=0.3067567046E+02,
c    .a53=0.4942041527E+01)
c    vbendb mar08/95    
c     parameter( 
c    .g11=-.4677295200E+02,g12=0.3587873970E+01,g13=-.7051720560E-01,
c    .g21=0.9549632701E-02,g22=0.9971623818E-02,g23=-.4040392363E-02,
c    .g24=0.9045746595E-03,g31=-.2746842691E+00,g32=-.1412945503E-01,
c    .g41=-.1337907240E-01,g42=-.6200741002E-03,g43=0.1257251291E+00,
c    .g44=-.4195584422E-01,g51=0.5287171694E+02,g52=0.4490366264E+01,
c    .g53=0.1735452866E+01)
c
c cbenda terms jun21/95                                                         
c-sngl;      parameter( 
c-sngl;     .c11=0.1860299931E+04,c12=-.6134458037E+03,c13=0.7337207161E+02,
c-sngl;     .c14=-.2676717625E+04,c15=0.1344099415E+04,c21=0.1538913137E+03,
c-sngl;     .c22=0.4348007369E+02,c23=0.1719720677E+03,c24=0.2115963042E+03,
c-sngl;     .c31=-.7026089414E+02,c32=-.1300938992E+03,c41=0.1310273564E+01,
c-sngl;     .c42=-.6175149574E+00,c43=-.2679089358E+02,c44=0.5577477171E+01,
c-sngl;     .c51=-.3543353539E+04,c52=-.3740709591E+03,c53=0.7979303144E+02,
c-sngl;     .c61=-.1104230585E+04,c62=0.4603572025E+04,c63=-.5593496634E+04,
c-sngl;     .c71=-.1069406434E+02,c72=0.1021807153E+01,c73=0.6669828341E-01,
c-sngl;     .c74=0.4168542348E+02,c75=0.1751608567E+02,c81=0.9486883238E+02,
c-sngl;     .c82=-.1519334221E+02,c83=0.4024697252E+04,c84=-.2225159395E+02) !jun21
      parameter( 
     .c11=0.1860299931d+04,c12=-.6134458037d+03,c13=0.7337207161d+02,
     .c14=-.2676717625d+04,c15=0.1344099415d+04,c21=0.1538913137d+03,
     .c22=0.4348007369d+02,c23=0.1719720677d+03,c24=0.2115963042d+03,
     .c31=-.7026089414d+02,c32=-.1300938992d+03,c41=0.1310273564d+01,
     .c42=-.6175149574d+00,c43=-.2679089358d+02,c44=0.5577477171d+01,
     .c51=-.3543353539d+04,c52=-.3740709591d+03,c53=0.7979303144d+02,
     .c61=-.1104230585d+04,c62=0.4603572025d+04,c63=-.5593496634d+04,
     .c71=-.1069406434d+02,c72=0.1021807153d+01,c73=0.6669828341d-01,
     .c74=0.4168542348d+02,c75=0.1751608567d+02,c81=0.9486883238d+02,
     .c82=-.1519334221d+02,c83=0.4024697252d+04,c84=-.2225159395d+02) !jun21
c           (last 4 parameters renumbered)
c cbendb terms jun21/95                                                         
c-sngl;      parameter( 
c-sngl;     .d11=0.4203543357E+03,d12=-.4922474096E+02,d13=0.3362942544E+00,
c-sngl;     .d14=-.3827423082E+03,d15=0.1746726001E+03,d21=0.1699995737E-01,
c-sngl;     .d22=0.1513036778E-01,d23=0.2659119354E-01,d24=-.5760387483E-02,
c-sngl;     .d31=0.1020622621E+02,d32=0.1050536271E-01,d41=0.6836172780E+00,
c-sngl;     .d42=-.1627858240E+00,d43=-.6925485045E+01,d44=0.1632567385E+01,
c-sngl;     .d51=0.1083595009E+04,d52=0.4641431791E+01,d53=-.8233144461E+00,
c-sngl;     .d61=-.6157225942E+02,d62=0.3094361471E+03,d63=-.3299631143E+03,
c-sngl;     .d71=0.8866227120E+01,d72=-.1382126854E+01,d73=0.7620770145E-01,
c-sngl;     .d74=-.5145757859E+02,d75=0.2046097265E+01,d81=0.2540775558E+01,
c-sngl;     .d82=-.4889246569E+00,d83=-.1127439280E+04,d84=-.2269932295E+01) !jun21
      parameter( 
     .d11=0.4203543357d+03,d12=-.4922474096d+02,d13=0.3362942544d+00,
     .d14=-.3827423082d+03,d15=0.1746726001d+03,d21=0.1699995737d-01,
     .d22=0.1513036778d-01,d23=0.2659119354d-01,d24=-.5760387483d-02,
     .d31=0.1020622621d+02,d32=0.1050536271d-01,d41=0.6836172780d+00,
     .d42=-.1627858240d+00,d43=-.6925485045d+01,d44=0.1632567385d+01,
     .d51=0.1083595009d+04,d52=0.4641431791d+01,d53=-.8233144461d+00,
     .d61=-.6157225942d+02,d62=0.3094361471d+03,d63=-.3299631143d+03,
     .d71=0.8866227120d+01,d72=-.1382126854d+01,d73=0.7620770145d-01,
     .d74=-.5145757859d+02,d75=0.2046097265d+01,d81=0.2540775558d+01,
     .d82=-.4889246569d+00,d83=-.1127439280d+04,d84=-.2269932295d+01) !jun21
c           (last 4 parameters renumbered)
c 
c    cbenda mar08/95   
c     parameter( 
c    .c11=0.1983833377E+04,c12=-.7161674985E+03,c13=0.9480354622E+02,
c    .c14=-.2860199829E+04,c15=0.1424701828E+04,c21=0.1544592401E+03,
c    .c22=0.3684848293E+02,c23=0.1717816399E+03,c24=0.2060270139E+03,
c    .c31=-.1215015279E+03,c32=-.1370574974E+03,c41=0.8903858200E+00,
c    .c42=-.6496543267E+00,c43=-.2367464822E+02,c44=0.5098769966E+01,
c    .c51=0.6728307258E+03,c52=-.3825078986E+03,c53=0.8128185587E+02,
c    .c61=-.1141379097E+04,c62=0.4771127549E+04,c63=-.5713570141E+04,
c    .c71=-.7394045711E+01,c72=0.9888694698E+00,c73=0.7493009017E-01,
c    .c74=0.9590644736E+02,c75=0.1925781563E+02,c81=0.1001295468E+03,
c    .c82=-.1533150643E+02,c83=-.1850297943E+03,c84=-.2420454635E+02)
c     (last 4 parameters renumbered)
c    cbendb terms mar08/95    
c     parameter( 
c    .d11=0.4441478933E+03,d12=-.5345816836E+02,d13=0.1443665553E+01,
c    .d14=-.4155134339E+03,d15=0.1912235377E+03,d21=-.3697316400E-02,
c    .d22=0.1688497284E-01,d23=0.2330055461E-01,d24=-.4670024282E-02,
c    .d31=0.6492099426E+01,d32=0.1597054504E-01,d41=0.6701003237E+00,
c    .d42=-.1868909419E+00,d43=-.6709596498E+01,d44=0.1713120242E+01,
c    .d51=-.5022064388E+04,d52=0.2322997491E+01,d53=-.2710031511E+00,
c    .d61=-.7265968242E+02,d62=0.3593588680E+03,d63=-.3933599161E+03,
c    .d71=0.9312252325E+01,d72=-.1380396384E+01,d73=0.7373026488E-01,
c    .d74=-.5070954743E+02,d75=0.7931122062E+00,d81=0.3599502276E+01,
c    .d82=-.5609868931E+00,d83=0.4980220561E+04,d84=-.1087179983E+01)
c     (last 4 parameters renumbered)
c
c Vbenda values from fit634b (used for surf706) 
c     parameter( 
c    .a11=-.2918252280E+03,a12=0.2164569141E+02,a13=-.4005497543E+00,
c    .a21=-.2890774947E+01,a22=0.1032032542E+02,a23=0.2681748056E+02,
c    .a24=0.2633751055E+02,a31=0.6180641351E+01,a32=0.5037667041E+01,
c    .a41=-.1125570079E+00,a42=-.3176529304E-02,a43=0.9068915355E+00,
c    .a44=-.7228418516E+00,a51=0.2785898232E+03,a52=0.4764442446E+02,
c    .a53=0.8621545662E+01)
c   Vbendb values from fit634b (used for surf706)
c     parameter( 
c    .g11=-.4241912309E+02,g12=0.2951365281E+01,g13=-.4840201514E-01,
c    .g21=-.1159168549E-03,g22=0.8688003567E-02,g23=-.3486923900E-02,
c    .g24=0.8312212839E-03,g31=0.5621810473E-01,g32=-.9776930747E-02,
c    .g41=-.1178456251E-01,g42=0.3491086729E-02,g43=0.7430516993E-01,
c    .g44=-.9643636957E-01,g51=0.4735533782E+02,g52=0.3001038808E+01,
c    .g53=0.1896630453E+01)
c
c parameters required for Cbend calculations:
c
c-- Cbenda values from fit706d.out 
c     parameter( 
c    .c11=0.7107064647E+04,c12=-.3728421267E+04,c13=0.1757654624E+03,
c    .c14=-.9725998132E+04,c15=0.4665086074E+04,c21=-.4435165986E+02,
c    .c22=0.1604477309E+03,c23=0.5805142925E+03,c24=0.6892349445E+03,
c    .c31=0.6581730442E+03,c32=0.6078389739E+02,c41=0.2885182566E+01,
c    .c42=-.1728169916E+01,c43=-.1119535503E+03,c44=0.4052536250E+02,
c    .c51=0.2540505673E+03,c52=-.5762083627E+03,c53=0.1295901320E+03,
c    .c61=-.2131706075E+04,c62=0.9084452020E+04,c63=-.1138253963E+05,
c    .c71=-.3964298833E+02,c72=-.5019979693E+01,c73=0.2906541488E+00,
c    .c74=0.1212943686E+04,c75=0.4140463415E+02,c81=0.1752855549E+03,
c    .c82=-.2496320107E+02,c83=0.3765413052E+03,c84=-.5480488130E+02)
c-- Cbendb values from fit706d.out 
c     parameter( 
c    .d11=0.1917166552E+04,d12=-.6542563392E+03,d13=0.6793758367E+02,
c    .d14=-.1694968847E+04,d15=0.6866649703E+03,d21=-.2137567948E+00,
c    .d22=0.4975938228E-01,d23=0.9364998295E-01,d24=-.2444320779E-01,
c    .d31=-.2863126914E+02,d32=0.5443219625E-01,d41=0.9673956120E+00,
c    .d42=-.1160159706E+01,d43=-.2424199759E+02,d44=0.8569424490E+01,
c    .d51=-.6517635862E+04,d52=0.1518098147E+03,d53=-.2706514366E+02,
c    .d61=0.4308392956E+02,d62=-.1234851732E+03,d63=0.2320626055E+03,
c    .d71=0.1049541418E+02,d72=-.2424169341E+01,d73=0.1745646946E+00,
c    .d74=0.2603615561E+02,d75=0.1345799970E+02,d81=-.6653710513E+01,
c    .d82=-.2576854447E+00,d83=0.6172608425E+04,d84=-.1328142473E+02)
c end of parameters
c
c   feb27/91 new c51 = c51+c83;   new d51 = d51+d83
      cx1 = c51 + c83
      dx1 = d51 + d83
      r1 = rpass(1)
      r2 = rpass(2)
      r3 = rpass(3)
      t1 = r1*r1 -r2*r2 -r3*r3 
      t2 = r2*r2 -r3*r3 -r1*r1 
      t3 = r3*r3 -r1*r1 -r2*r2 
c  calculate the cosines of the three internal angles:
      c1 = t1 / (-2.d0*r2*r3)
      c2 = t2 / (-2.d0*r3*r1)
      c3 = t3 / (-2.d0*r1*r2)
      sum = c1 + c2 + c3
      B1a = 1.d0 - sum         
      c1cube = c1*c1*c1
      c2cube = c2*c2*c2
      c3cube = c3*c3*c3
      cos3t1 = 4.d0*c1cube - 3.d0*c1
      cos3t2 = 4.d0*c2cube - 3.d0*c2
      cos3t3 = 4.d0*c3cube - 3.d0*c3
      sumb   = cos3t1 + cos3t2 + cos3t3 
      B1b    = 1.d0 - ( z58*sumb + z38*sum )
c  calculate derivatives if desired
         dc1dr1 = -r1/(r2*r3)
         dc2dr2 = -r2/(r1*r3)
         dc3dr3 = -r3/(r1*r2)
         dc1dr2 = ( t1/(r2*r2) + 2.d0 )/(2.d0*r3)
         dc1dr3 = ( t1/(r3*r3) + 2.d0 )/(2.d0*r2)
         dc2dr1 = ( t2/(r1*r1) + 2.d0 )/(2.d0*r3)
         dc2dr3 = ( t2/(r3*r3) + 2.d0 )/(2.d0*r1)
         dc3dr1 = ( t3/(r1*r1) + 2.d0 )/(2.d0*r2)
         dc3dr2 = ( t3/(r2*r2) + 2.d0 )/(2.d0*r1)
         db1a(1) = -1.d0*( dc1dr1 + dc2dr1 + dc3dr1 )
         db1a(2) = -1.d0*( dc1dr2 + dc2dr2 + dc3dr2 )
         db1a(3) = -1.d0*( dc1dr3 + dc2dr3 + dc3dr3 )
         d1 = 12.d0*c1*c1 - 3.d0
         d2 = 12.d0*c2*c2 - 3.d0
         d3 = 12.d0*c3*c3 - 3.d0
         db1b(1) = -z58*( d1*dc1dr1 + d2*dc2dr1 + d3*dc3dr1 )
     .             -z38*(    dc1dr1 +    dc2dr1 +    dc3dr1 )
         db1b(2) = -z58*( d1*dc1dr2 + d2*dc2dr2 + d3*dc3dr2 )
     .             -z38*(    dc1dr2 +    dc2dr2 +    dc3dr2 )
         db1b(3) = -z58*( d1*dc1dr3 + d2*dc2dr3 + d3*dc3dr3 )
     .             -z38*(    dc1dr3 +    dc2dr3 +    dc3dr3 )
cd    if(ipr.gt.1)then
cd        write(7,*)
cd        write(7,*) ' from subr.vbcb -----------------------------'
cd        write(7,6000) rpass
cd        write(7,6010) b1a,b1b
cd        write(7,6020) db1a,db1b
cd        if(ipr.gt.2)then
cd           write(7,7000) dc1dr1,dc1dr2,dc1dr3,
cd   .                     dc2dr1,dc2dr2,dc2dr3,
cd   .                     dc3dr1,dc3dr2,dc3dr3
cd           write(7,7010) t1,t2,t3,c1,c2,c3,d1,d2,d3
cd        end if
cd    end if
c
c  calculate the quantities used by both Vbend and Cbend:
c
      R    = r1 + r2 + r3
      Rsq  = R*R
      B2   = 1.d0/r1 + 1.d0/r2 + 1.d0/r3
      B3   = (r2-r1)*(r2-r1) + (r3-r2)*(r3-r2) + (r1-r3)*(r1-r3)
      eps2 = 1.d-12
      B3b  = dsqrt( B3 + eps2 )
	 dB2(1) = -1.d0/(r1*r1)
	 dB2(2) = -1.d0/(r2*r2)
	 dB2(3) = -1.d0/(r3*r3)
	 dB3(1) = 4.d0*r1 -2.d0*r2 -2.d0*r3
	 dB3(2) = 4.d0*r2 -2.d0*r3 -2.d0*r1
	 dB3(3) = 4.d0*r3 -2.d0*r1 -2.d0*r2
	 dB3(4) = 0.5d0*dB3(1)/B3b
	 dB3(5) = 0.5d0*dB3(2)/B3b
	 dB3(6) = 0.5d0*dB3(3)/B3b
         exp1   = dexp(-beta1*R)
         exp2   = dexp(-beta2*Rsq)
         exp7   = dexp(-beta3*R)
	 dexp1  = -beta1*exp1
	 dexp2  = -2.d0*beta2*R*exp2
	 dexp7  = -beta3*exp7
c   do the vbnda calculations:
         B1    = B1a
         B12   = B1*B1
         B13   = B12*B1
         B14   = B13*B1
         B15   = B14*B1
         asum  = a11 +  a12*R +  a13*Rsq 
         bsum  = a21*B12 +  a22*B13 +  a23*B14 +  a24*B15
	 csum  = a31*B1*exp1 +  a32*B12*exp2
         dsum1 = a41*exp1 +  a42*exp2
         dsum2 = a43*exp1 +  a44*exp2
         fsum  = a51 +  a52*R +  a53*Rsq 
         vb(1,1) = B1*asum*exp1 
         vb(1,2) = bsum * exp2
         vb(1,3) = B2*csum
         vb(1,4) = B1*B3*dsum1 + B1*B3b*dsum2
         vb(1,5) = B1* fsum  *exp7 
         Vbnd(1) = vb(1,1)+vb(1,2)+vb(1,3)+vb(1,4)+vb(1,5)
c
	 dasum   = a12+2.d0*a13*R
	 dbsum   =  2.d0* a21*B1  + 3.d0* a22*B12
     .            + 4.d0* a23*B13 + 5.d0* a24*B14
	 ddsum1  =  a41*dexp1 +  a42*dexp2
	 ddsum2  =  a43*dexp1 +  a44*dexp2
	 dfsum   =  a52 +2.d0*a53*R
         dVb1(1) = dB1a(1)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1
         dVb1(2) = dB1a(2)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1
         dVb1(3) = dB1a(3)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1
c 
         dVb2(1) = dbsum*dB1a(1)*exp2 + bsum*dexp2
         dVb2(2) = dbsum*dB1a(2)*exp2 + bsum*dexp2
         dVb2(3) = dbsum*dB1a(3)*exp2 + bsum*dexp2
c
c calculate the Vb3 derivatives:
         dVb3(1)=dB2(1)*csum + B2*(  a31*dB1a(1)*exp1 +  a31*B1*dexp1
     .                  +2.d0* a32*B1*dB1a(1)*exp2 +  a32*B12*dexp2)
         dVb3(2)=dB2(2)*csum + B2*(  a31*dB1a(2)*exp1 +  a31*B1*dexp1
     .                  +2.d0* a32*B1*dB1a(2)*exp2 +  a32*B12*dexp2)
         dVb3(3)=dB2(3)*csum + B2*(  a31*dB1a(3)*exp1 +  a31*B1*dexp1
     .                  +2.d0* a32*B1*dB1a(3)*exp2 +  a32*B12*dexp2)
c
c calculate the Vb4 derivatives (may 27/90):
         dVb4(1)= dB1a(1)*B3 *dsum1+  B1*dB3(1)*dsum1 + B1*B3 *ddsum1
     .          + dB1a(1)*B3b*dsum2 + B1*dB3(4)*dsum2 + B1*B3b*ddsum2
         dVb4(2)= dB1a(2)*B3 *dsum1+  B1*dB3(2)*dsum1 + B1*B3 *ddsum1
     .          + dB1a(2)*B3b*dsum2 + B1*dB3(5)*dsum2 + B1*B3b*ddsum2
         dVb4(3)= dB1a(3)*B3 *dsum1+  B1*dB3(3)*dsum1 + B1*B3 *ddsum1
     .          + dB1a(3)*B3b*dsum2 + B1*dB3(6)*dsum2 + B1*B3b*ddsum2
c
c calculate the Vb5 derivatives:
	 dVb5(1) = dB1a(1)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7
	 dVb5(2) = dB1a(2)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7
	 dVb5(3) = dB1a(3)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7
c
c calculate the overall derivatives:
         dVbnda(1) = dvb1(1)+dvb2(1)+dvb3(1)+dvb4(1)+dvb5(1)
         dVbnda(2) = dvb1(2)+dvb2(2)+dvb3(2)+dvb4(2)+dvb5(2)
         dVbnda(3) = dvb1(3)+dvb2(3)+dvb3(3)+dvb4(3)+dvb5(3)
	 Vbnda = Vbnd(1)
c
c---- now do the vbendb calculations --------
         B1    = B1b
         B12   = B1*B1
         B13   = B12*B1
         B14   = B13*B1
         B15   = B14*B1
         asum  = g11 +  g12*R +  g13*Rsq 
         bsum  = g21*B12 +  g22*B13 +  g23*B14 +  g24*B15
	 csum  = g31*B1*exp1 +  g32*B12*exp2
         dsum1 = g41*exp1 +  g42*exp2
         dsum2 = g43*exp1 +  g44*exp2
         fsum  = g51 +  g52*R +  g53*Rsq 
         vb(2,1) = B1*asum*exp1 
         vb(2,2) = bsum * exp2
         vb(2,3) = B2*csum
         vb(2,4) = B1*B3*dsum1 + B1*B3b*dsum2
         vb(2,5) = B1* fsum  *exp7 
         Vbnd(2) = vb(2,1)+vb(2,2)+vb(2,3)+vb(2,4)+vb(2,5)
c
	 dasum   = g12+2.d0*g13*R
	 dbsum   =  2.d0* g21*B1  + 3.d0* g22*B12
     .             +4.d0* g23*B13 + 5.d0* g24*B14
	 ddsum1  =  g41*dexp1 +  g42*dexp2
	 ddsum2  =  g43*dexp1 +  g44*dexp2
	 dfsum   =  g52 +2.d0*g53*R
         dVb1(1) = dB1b(1)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1
         dVb1(2) = dB1b(2)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1
         dVb1(3) = dB1b(3)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1
c 
         dVb2(1) = dbsum*dB1b(1)*exp2 + bsum*dexp2
         dVb2(2) = dbsum*dB1b(2)*exp2 + bsum*dexp2
         dVb2(3) = dbsum*dB1b(3)*exp2 + bsum*dexp2
c
c calculate the Vb3 derivatives:
         dVb3(1)=dB2(1)*csum + B2*(  g31*dB1b(1)*exp1 +  g31*B1*dexp1
     .                  +2.d0* g32*B1*dB1b(1)*exp2 +  g32*B12*dexp2)
         dVb3(2)=dB2(2)*csum + B2*(  g31*dB1b(2)*exp1 +  g31*B1*dexp1
     .                  +2.d0* g32*B1*dB1b(2)*exp2 +  g32*B12*dexp2)
         dVb3(3)=dB2(3)*csum + B2*(  g31*dB1b(3)*exp1 +  g31*B1*dexp1
     .                  +2.d0* g32*B1*dB1b(3)*exp2 +  g32*B12*dexp2)
c
c calculate the Vb4 derivatives (may 27/90):
         dVb4(1)= dB1b(1)*B3 *dsum1+  B1*dB3(1)*dsum1 + B1*B3 *ddsum1
     .          + dB1b(1)*B3b*dsum2 + B1*dB3(4)*dsum2 + B1*B3b*ddsum2
         dVb4(2)= dB1b(2)*B3 *dsum1+  B1*dB3(2)*dsum1 + B1*B3 *ddsum1
     .          + dB1b(2)*B3b*dsum2 + B1*dB3(5)*dsum2 + B1*B3b*ddsum2
         dVb4(3)= dB1b(3)*B3 *dsum1+  B1*dB3(3)*dsum1 + B1*B3 *ddsum1
     .          + dB1b(3)*B3b*dsum2 + B1*dB3(6)*dsum2 + B1*B3b*ddsum2
c
c calculate the Vb5 derivatives:
	 dVb5(1) = dB1b(1)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7
	 dVb5(2) = dB1b(2)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7
	 dVb5(3) = dB1b(3)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7
c
c calculate the overall derivatives:
         dVbndb(1) = dvb1(1)+dvb2(1)+dvb3(1)+dvb4(1)+dvb5(1)
         dVbndb(2) = dvb1(2)+dvb2(2)+dvb3(2)+dvb4(2)+dvb5(2)
         dVbndb(3) = dvb1(3)+dvb2(3)+dvb3(3)+dvb4(3)+dvb5(3)
	 Vbndb = Vbnd(2)
c
c  now calculate the Cbend correction terms:
c  Cbnda uses the B1a formula and Cbndb uses the B1b formula
      if(icompc.eq.0) return
      sumT  = T(1) + T(2) + T(3)
      Rcu   = Rsq*R
      P     = r1*r2*r3
      Psq   = P*P
      Pcu   = Psq*P
      dP(1) = r2*r3
      dP(2) = r3*r1
      dP(3) = r1*r2
      Cbnds(1)  = 0.d0
      Cbnds(2)  = 0.d0
      dCbnda(1) = 0.d0
      dCbnda(2) = 0.d0
      dCbnda(3) = 0.d0
      dCbndb(1) = 0.d0
      dCbndb(2) = 0.d0
      dCbndb(3) = 0.d0
c  calculate exponentials and derivatives (common to cbnda and cbndb):
            exp7  = dexp( -beta2*Rcu ) 
	    dexp7 = -3.d0*beta2*Rsq*exp7
c ----- calculate the Cbnda correction term ----------------
       B1    = B1a
       B12   = B1*B1
       B13   = B12*B1
       B14   = B13*B1
       B15   = B14*B1
       asum  = c11 + c12*R + c13*Rsq + c14/R + c15/Rsq 
       bsum  = c21*B12 + c22*B13 + c23*B14 + c24*B15
       csum  = c31*B1 *exp1  + c32*B12*exp2  
       dsum1 = c41*exp1 + c42*exp2
       dsum2 = c43*exp1 + c44*exp2
       fsum  = cx1 + c52*R + c53*Rsq 
       gsum  = c61 + c62/R + c63/Rsq 
       aasum = c71 + c72*P + c73*Psq + c74/P + c75/Psq 
       ffsum =       c81*P + c82*Psq + c84/Psq 
       dasum = c12 + 2.d0*c13*R -c14/Rsq -2.d0*c15/Rcu
       dbsum = 2.d0*c21*B1 +3.d0*c22*B12 +4.d0*c23*B13 +5.d0*c24*B14
       ddsum1= c41*dexp1 + c42*dexp2
       ddsum2= c43*dexp1 + c44*dexp2
       dfsum = c52 + 2.d0*c53*R
       dgsum = -c62/Rsq - 2.d0*c63/Rcu
       daasum = c72 +2.d0*c73*P -c74/Psq -2.d0*c75/Pcu
       dffsum = c81 + 2.d0*c82*P -2.d0*c84/Pcu
         cb(1,1)  = B1*asum * exp1  /P
         cb(1,2)  = bsum * exp2  
         cb(1,3)  = B2*csum 
         cb(1,4)  = B1*B3*dsum1 + B1*B3b*dsum2
         cb(1,5)  = B1 * fsum * exp7 /P
         cb(1,6)  = B1 * gsum * exp7  
         cb(1,7)  = B1 * aasum * exp2 
         cb(1,8)  = B1 * ffsum * exp7 
         Cbnds(1) = cb(1,1)+cb(1,2)+cb(1,3)+cb(1,4)+cb(1,5)
     .                     +cb(1,6)+cb(1,7)+cb(1,8)
c       calculate the derivatives:
	 dCb1(1) = dB1a(1)*asum*exp1/P + B1*dasum*exp1/P
     .                +B1*asum*dexp1/P - B1*asum*exp1*dP(1)/Psq
	 dCb1(2) = dB1a(2)*asum*exp1/P + B1*dasum*exp1/P
     .                +B1*asum*dexp1/P - B1*asum*exp1*dP(2)/Psq
	 dCb1(3) = dB1a(3)*asum*exp1/P + B1*dasum*exp1/P
     .                +B1*asum*dexp1/P - B1*asum*exp1*dP(3)/Psq
c
	 dCb2(1) = dbsum*dB1a(1)*exp2 + bsum*dexp2
	 dCb2(2) = dbsum*dB1a(2)*exp2 + bsum*dexp2
	 dCb2(3) = dbsum*dB1a(3)*exp2 + bsum*dexp2
c
	 dCb3(1) = dB2(1)*csum + B2*( c31*dB1a(1)*exp1 +c31*B1*dexp1
     .            +c32 *2.d0*B1*dB1a(1)*exp2 + c32 *B12*dexp2 )
	 dCb3(2) = dB2(2)*csum + B2*( c31 *dB1a(2)*exp1 +c31 *B1*dexp1
     .            +c32 *2.d0*B1*dB1a(2)*exp2 + c32 *B12*dexp2 )
	 dCb3(3) = dB2(3)*csum + B2*( c31 *dB1a(3)*exp1 +c31 *B1*dexp1
     .            +c32 *2.d0*B1*dB1a(3)*exp2 + c32 *B12*dexp2 )
c          dCb4 equations may 27/90:
	 dCb4(1) = dB1a(1)*B3 *dsum1 + B1*dB3(1)*dsum1 + B1*B3 *ddsum1
     .           + dB1a(1)*B3b*dsum2 + B1*dB3(4)*dsum2 + B1*B3b*ddsum2
	 dCb4(2) = dB1a(2)*B3 *dsum1 + B1*dB3(2)*dsum1 + B1*B3 *ddsum1
     .           + dB1a(2)*B3b*dsum2 + B1*dB3(5)*dsum2 + B1*B3b*ddsum2
	 dCb4(3) = dB1a(3)*B3 *dsum1 + B1*dB3(3)*dsum1 + B1*B3 *ddsum1
     .           + dB1a(3)*B3b*dsum2 + B1*dB3(6)*dsum2 + B1*B3b*ddsum2
c          dCb5 equations may 27/90: (corrected jun01)
	 dCb5(1) = dB1a(1)*fsum*exp7/P + B1*dfsum*exp7/P
     .               + B1*fsum*dexp7/P - B1*fsum*exp7*dP(1)/Psq
	 dCb5(2) = dB1a(2)*fsum*exp7/P + B1*dfsum*exp7/P
     .               + B1*fsum*dexp7/P - B1*Fsum*exp7*dP(2)/Psq
	 dCb5(3) = dB1a(3)*fsum*exp7/P + B1*dfsum*exp7/P
     .               + B1*fsum*dexp7/P - B1*Fsum*exp7*dP(3)/Psq
c
	 dCb6(1) = dB1a(1)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7
	 dCb6(2) = dB1a(2)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7
	 dCb6(3) = dB1a(3)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7
c
	 dCb7(1) = dB1a(1)*aasum*exp2 + B1*daasum*dP(1)*exp2
     .                               + B1* aasum*dexp2
	 dCb7(2) = dB1a(2)*aasum*exp2 + B1*daasum*dP(2)*exp2
     .                               + B1* aasum*dexp2
	 dCb7(3) = dB1a(3)*aasum*exp2 + B1*daasum*dP(3)*exp2
     .                               + B1* aasum*dexp2
	 dCb8(1) = dB1a(1)*ffsum*exp7  +B1*dffsum*dP(1)*exp7 
     .                                +B1* ffsum*dexp7 
	 dCb8(2) = dB1a(2)*ffsum*exp7  +B1*dffsum*dP(2)*exp7 
     .                                +B1* ffsum*dexp7 
	 dCb8(3) = dB1a(3)*ffsum*exp7  +B1*dffsum*dP(3)*exp7 
     .                                +B1* ffsum*dexp7 
           dCbnds(1,1) = dCb1(1) +dCb2(1) +dCb3(1) +dCb4(1) +dCb5(1)
     .                  +dCb6(1) +dCb7(1) +dCb8(1) 
           dCbnds(1,2) = dCb1(2) +dCb2(2) +dCb3(2) +dCb4(2) +dCb5(2)
     .                  +dCb6(2) +dCb7(2) +dCb8(2) 
           dCbnds(1,3) = dCb1(3) +dCb2(3) +dCb3(3) +dCb4(3) +dCb5(3)
     .                  +dCb6(3) +dCb7(3) +dCb8(3) 
c
c ----- calculate the Cbndb correction term ----------------
       B1    = B1b
       B12   = B1*B1
       B13   = B12*B1
       B14   = B13*B1
       B15   = B14*B1
       asum  = d11 + d12*R + d13*Rsq + d14/R + d15/Rsq 
       bsum  = d21*B12 + d22*B13 + d23*B14 + d24*B15
       csum  = d31*B1 *exp1  + d32*B12*exp2  
       dsum1 = d41*exp1 + d42*exp2
       dsum2 = d43*exp1 + d44*exp2
       fsum  = dx1 + d52*R + d53*Rsq 
       gsum  = d61 + d62/R + d63/Rsq 
       aasum = d71 + d72*P + d73*Psq + d74/P + d75/Psq 
       ffsum =       d81*P + d82*Psq + d84/Psq 
       dasum = d12 + 2.d0*d13*R -d14/Rsq -2.d0*d15/Rcu
       dbsum = 2.d0*d21*B1 +3.d0*d22*B12 +4.d0*d23*B13 +5.d0*d24*B14
       ddsum1= d41*dexp1 + d42*dexp2
       ddsum2= d43*dexp1 + d44*dexp2
       dfsum = d52 + 2.d0*d53*R
       dgsum = -d62/Rsq - 2.d0*d63/Rcu
       daasum = d72 +2.d0*d73*P -d74/Psq -2.d0*d75/Pcu
       dffsum = d81 + 2.d0*d82*P -2.d0*d84/Pcu
         cb(2,1)  = B1*asum * exp1  /P
         cb(2,2)  = bsum * exp2  
         cb(2,3)  = B2*csum 
         cb(2,4)  = B1*B3*dsum1 + B1*B3b*dsum2
         cb(2,5)  = B1 * fsum * exp7 /P
         cb(2,6)  = B1 * gsum * exp7  
         cb(2,7)  = B1 * aasum * exp2 
         cb(2,8)  = B1 * ffsum * exp7  
         Cbnds(2) = cb(2,1)+cb(2,2)+cb(2,3)+cb(2,4)+cb(2,5)
     .                     +cb(2,6)+cb(2,7)+cb(2,8)
c       calculate the derivatives:
	 dCb1(1) = dB1b(1)*asum*exp1/P + B1*dasum*exp1/P
     .                +B1*asum*dexp1/P - B1*asum*exp1*dP(1)/Psq
	 dCb1(2) = dB1b(2)*asum*exp1/P + B1*dasum*exp1/P
     .                +B1*asum*dexp1/P - B1*asum*exp1*dP(2)/Psq
	 dCb1(3) = dB1b(3)*asum*exp1/P + B1*dasum*exp1/P
     .                +B1*asum*dexp1/P - B1*asum*exp1*dP(3)/Psq
c
	 dCb2(1) = dbsum*dB1b(1)*exp2 + bsum*dexp2
	 dCb2(2) = dbsum*dB1b(2)*exp2 + bsum*dexp2
	 dCb2(3) = dbsum*dB1b(3)*exp2 + bsum*dexp2
c
	 dCb3(1)= dB2(1)*csum + B2*( d31 *dB1b(1)*exp1 +d31 *B1*dexp1
     .           +d32 *2.d0*B1*dB1b(1)*exp2 + d32 *B12*dexp2 )
	 dCb3(2)= dB2(2)*csum + B2*( d31 *dB1b(2)*exp1 +d31 *B1*dexp1
     .           +d32 *2.d0*B1*dB1b(2)*exp2 + d32 *B12*dexp2 )
	 dCb3(3)= dB2(3)*csum + B2*( d31 *dB1b(3)*exp1 +d31 *B1*dexp1
     .           +d32 *2.d0*B1*dB1b(3)*exp2 + d32 *B12*dexp2 )
c          dCb4 equations may 27/90:
	 dCb4(1) = dB1b(1)*B3 *dsum1 + B1*dB3(1)*dsum1 + B1*B3 *ddsum1
     .           + dB1b(1)*B3b*dsum2 + B1*dB3(4)*dsum2 + B1*B3b*ddsum2
	 dCb4(2) = dB1b(2)*B3 *dsum1 + B1*dB3(2)*dsum1 + B1*B3 *ddsum1
     .           + dB1b(2)*B3b*dsum2 + B1*dB3(5)*dsum2 + B1*B3b*ddsum2
	 dCb4(3) = dB1b(3)*B3 *dsum1 + B1*dB3(3)*dsum1 + B1*B3 *ddsum1
     .           + dB1b(3)*B3b*dsum2 + B1*dB3(6)*dsum2 + B1*B3b*ddsum2
c          dCb5 equations may 27/90:
	 dCb5(1) = dB1b(1)*fsum*exp7/P + B1*dfsum*exp7/P
     .               + B1*fsum*dexp7/P - B1*fsum*exp7*dP(1)/Psq
	 dCb5(2) = dB1b(2)*fsum*exp7/P + B1*dfsum*exp7/P
     .               + B1*fsum*dexp7/P - B1*fsum*exp7*dP(2)/Psq
	 dCb5(3) = dB1b(3)*fsum*exp7/P + B1*dfsum*exp7/P
     .               + B1*fsum*dexp7/P - B1*fsum*exp7*dP(3)/Psq
c
	 dCb6(1) = dB1b(1)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7
	 dCb6(2) = dB1b(2)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7
	 dCb6(3) = dB1b(3)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7
c
	 dCb7(1) = dB1b(1)*aasum*exp2 + B1*daasum*dP(1)*exp2
     .                               + B1* aasum*dexp2
	 dCb7(2) = dB1b(2)*aasum*exp2 + B1*daasum*dP(2)*exp2
     .                               + B1* aasum*dexp2
	 dCb7(3) = dB1b(3)*aasum*exp2 + B1*daasum*dP(3)*exp2
     .                               + B1* aasum*dexp2
	 dCb8(1) = dB1b(1)*ffsum*exp7  +B1*dffsum*dP(1)*exp7 
     .                                +B1* ffsum*dexp7 
	 dCb8(2) = dB1b(2)*ffsum*exp7  +B1*dffsum*dP(2)*exp7 
     .                                +B1* ffsum*dexp7 
	 dCb8(3) = dB1b(3)*ffsum*exp7  +B1*dffsum*dP(3)*exp7 
     .                                +B1* ffsum*dexp7 
           dCbnds(2,1) = dCb1(1) +dCb2(1) +dCb3(1) +dCb4(1) +dCb5(1)
     .                  +dCb6(1) +dCb7(1) +dCb8(1) 
           dCbnds(2,2) = dCb1(2) +dCb2(2) +dCb3(2) +dCb4(2) +dCb5(2)
     .                  +dCb6(2) +dCb7(2) +dCb8(2) 
           dCbnds(2,3) = dCb1(3) +dCb2(3) +dCb3(3) +dCb4(3) +dCb5(3)
     .                  +dCb6(3) +dCb7(3) +dCb8(3) 
c
c  calculate the total derivative from the pieces:
         dCbnda(1) = dT(1) * Cbnds(1)    + sumT * dCbnds(1,1)
         dCbndb(1) = dT(1) * Cbnds(2)    + sumT * dCbnds(2,1) 
         dCbnda(2) = dT(2) * Cbnds(1)    + sumT * dCbnds(1,2)
         dCbndb(2) = dT(2) * Cbnds(2)    + sumT * dCbnds(2,2) 
         dCbnda(3) = dT(3) * Cbnds(1)    + sumT * dCbnds(1,3)
         dCbndb(3) = dT(3) * Cbnds(2)    + sumT * dCbnds(2,3) 
      Cbnda = sumT*Cbnds(1)
      Cbndb = sumT*Cbnds(2)
c  some debugging print statements:
cd    if(ipr.gt.1)then
cd       write(7,7650) B2,B3
cd       write(7,7700) db2,db3
cd       write(7,*)
cd       write(7,*) 'Vbend values :'
cd       k=1
cd       write(7,7150) k,dvb1,dvb2,dvb3,dvb4,dvb5,dVbnda
cd       write(7,7155) (vb(1,i),i=1,5)
cd       k=2
cd       write(7,7150) k,dvb1,dvb2,dvb3,dvb4,dvb5,dVbndb
cd       write(7,7155) (vb(2,i),i=1,5)
cd       write(7,7100) Vbnda,Vbndb
cd       write(7,*)
cd       write(7,*) 'Cbend values:'
cd       write(7,*) 'c81,c81,c82=',c81,c81,c82
cd       write(7,*) 'c91,c83,c84=',c91,c83,c84
cd       write(7,*) 'd81,d81,d82=',d81,d81,d82
cd       write(7,*) 'd91,d83,d84=',d91,d83,d84
cd       k=1              
cd          write(7,2020) k,(cb(1,j),j=1,9)
cd          write(7,2033) sumT,Cbnds(1)
cd       k=2              
cd          write(7,2020) k,(cb(k,j),j=1,9)
cd          write(7,2033) sumT,Cbnds(2)
cd       k=1
cd       write(7,7152) k,dcb1,dcb2,dcb3,dcb4,dcb5,dcb6,
cd   .               dcb7,dcb8,dcb9,(dCbnds(1,i),i=1,3)
cd       k=2
cd       write(7,7152) k,dcb1,dcb2,dcb3,dcb4,dcb5,dcb6,
cd   .                 dcb7,dcb8,dcb9,(dCbnds(k,i),i=1,3)
cd       write(7,7200) dT,dCbnda,dCbndb
cd       write(7,*) 'VbCb  exit -------------------------'
cd       write(7,*)
cd    end if
      return
 2000 format(5x,'r1, r2, r3 = ',3(1x,f12.8))
 2010 format(5x,'B1a, B1b, B2, B3 = ',4(1x,g12.6))
 7100 format('    Vbnda = ',g12.6,'     Vbndb = ',g12.6)
 7650 format('    B2 =    ',g12.6,'       B3 =  ',g12.6)
 7700 format('    dB2 =   ',3(1x,g12.6),/,
     .       '    dB3 =   ',6(1x,g12.6))
 7150 format(' k=',i1,'     derivatives: ',/,
     .       ' dVb1  = ',3(1x,g16.10),/,
     .       ' dVb2  = ',3(1x,g16.10),/,
     .       ' dVb3  = ',3(1x,g16.10),/,
     .       ' dVb4  = ',3(1x,g16.10),/,
     .       ' dVb5  = ',3(1x,g16.10),/,
     .       ' dVbnd = ',3(1x,g16.10))
 7155 format('  Vb1 = ',g16.10,'  Vb2 = ',g16.10,/,
     .       '  Vb3 = ',g16.10,'  Vb4 = ',g16.10,/,
     .       '  Vb5 = ',g16.10)
 2020 format(5x,'k=',i1,' cb1,cb2,cb3= ',3(1x,g16.10),
     .     /,8x,' cb4,cb5,cb6= ',3(1x,g16.10),
     .     /,8x,' cb7,cb8,cb9= ',3(1x,g16.10))
 2033 format(5x,'sumT = ',g12.6 ,'  cbnds(k) = ',g12.6)
 7152 format('Subr.Cbend:   k=',i1,'     derivatives: ',/,
     .       '  dCb1  = ',3(1x,g18.12),/,
     .       '  dCb2  = ',3(1x,g18.12),/,
     .       '  dCb3  = ',3(1x,g18.12),/,
     .       '  dCb4  = ',3(1x,g18.12),/,
     .       '  dCb5  = ',3(1x,g18.12),/,
     .       '  dCb6  = ',3(1x,g18.12),/,
     .       '  dCb7  = ',3(1x,g18.12),/,
     .       '  dCb8  = ',3(1x,g18.12),/,
     .       '  dCb9  = ',3(1x,g18.12),/,
     .       '  dCbnd = ',3(1x,g18.12))
 7200 format('  dT    = ',3(1x,g18.12),/,
     .       '  dCbnda= ',3(1x,g18.12),/,
     .       '  dCbndb= ',3(1x,g18.12))
 6000 format(' from subr.b1ab:  r=',3(1x,f9.6))
 6010 format('    b1a = ',f16.10,'     b1b = ',f16.10)
 6020 format('    db1a: ',3(1x,e12.6),/,'    db1b: ',3(1x,e12.6))
 7000 format('    dc(j)/dr(i) derivatives: ',/,
     .    8x,3(1x,e12.6),/,8x,3(1x,e12.6),/,8x,3(1x,e12.6))
 7010 format('    ri2-rj2-rk2: ',3(1x,e12.6),/,
     .       '       cosines : ',3(1x,e12.6),/,
     .       '    12c**2 - 3 : ',3(1x,e12.6))
      end
c
      subroutine chgeom(r,ivalid)
c----------------------------------------------------------------------c
c apr21/95
c first, check that it is a valid H3 geometry, within tolerance "derr":
c this test corrected on nov23/93
c     generally, rlo + rmed > rhi
c     sometimes, rlo + rmed = rhi (linear geometry)
c     but if     rlo + rmed < rhi there's a problem
c
c error1 ... check that rlo+rmid>rhi
c istop1 ... [0] warning [1]stop calculations
c
c error2 ... check that all r(i)'s are > 0.2 bohr
c istop2 ... [0] warning [1]stop calculations
c----------------------------------------------------------------------c
      implicit double precision (a-h,o-z)
      parameter(derr = 1.d-5)
      parameter( istop1 = 1 , istop2 = 1 )
      dimension r(3)
      ivalid = 1
      rhi = dmax1(r(1),r(2))
      if(r(3).ge.rhi)then
         rmid = rhi
         rhi = r(3)
      else
         rmid = dmax1( r(3) , dmin1(r(1),r(2)) )
      endif
      rlo = dmin1(r(1),r(2),r(3))
cx    if(rlo+rmid.gt.rhi+derr)then  {old condition}
      if(rlo+rmid+derr.lt.rhi)then
cd       write(6,*) ' rlongest          = ',rhi+derr
cd       write(6,*) ' rmiddle,rshortest = ',rmid,rlo
         write(6,*) 'Warning: invalid geometry ',r
         if(istop1.gt.0) stop ' STOP -- geometry error '
c         ivalid = -1
      end if
      if(rlo.lt.0.2d0)then
cd       write(6,*) ' shortest distance: ',rlo
         write(6,*) 'Warning: invalid geometry ',r
         if(istop2.gt.0) stop ' STOP -- geometry error '
c         ivalid = -2
      end if
      return
      end
